Monte Carlo Quasi-heatbath by Approximate Inversion

When sampling the distribution P (~ φ) ∝ exp(−|A~ φ|2), a global heatbath normally proceeds by solving the linear system A~ φ = ~η, where ~η is a normal Gaussian vector, exactly. This paper shows how to preserve the distribution P (~ φ) while solving the linear system with arbitrarily low accuracy. Generalizations are presented. PACS numbers: 02.70.L, 02.50.N, 52.65.P, 11.15.H, 12.38.G In Monte Carlo simulations, it is frequently the case that one wants to sample a vector ~ φ from a distribution of the Gaussian type ∝ exp(−|A~ φ|2). Typically, ~ φ has many components, andA is a large, sparse matrix. In lattice field theory, ~ φ is the value of the continuum field ~ φ at regular grid points, and A is the discretized version of some differential operator A. Illustrative examples used in this paper are A = m + i~ p (free field) and A = m + ip/ (Dirac operator). The goal of the Monte Carlo simulation is to provide independent configurations of ~ φ at the least cost. The brute force approach consists in drawing successive random vectors ~η(k) from the normal Gaussian distribution exp(−|~η|2), and in solving A~ φ(k) = ~η(k). The solution of this linear system can be efficiently obtained with an iterative linear solver (Conjugate Gradient if A is Hermitian, BiCGStab otherwise). This approach can be called a global heatbath, because ~ φ(k+1) has no memory of ~ φ(k): the heatbath has touched all the components of ~ φ. To avoid a bias, the solver must be iterated to full convergence, which is often prohibitively expensive. One may try to limit the accuracy while maintaining the bias below statistical errors, but this requires a delicate compromise difficult to tune a priori. A notable example of this global heatbath approach is the stochastic evaluation of inverse Dirac matrix elements, where several hundred “noise vectors” ~η(k) are inverted to yield (A†A) ij ≈ 〈φiφ † j〉k. An abundant literature has been devoted to the optimization of this procedure [1, 2]. For the free field or the Dirac operator mentioned above, the number of iterations of the solver required to reach a given accuracy grows like the correlation length ξ ≡ 1/m. Thus the work per new, independent ~ φ is c ξz where z, the dynamical critical exponent, is 1. However the prefactor c is large. For this reason, local updates, where only one component of ~ φ is changed at a time, are often preferred. They usually provide an independent ~ φ after an amount of work c′ ξ2, but with a much smaller prefactor c′. 1 This paper presents an adaptation of the global heatbath, which allows for arbitrarily low accuracy in the solution of A~ φ = ~η, thus reducing the prefactor c, while maintaining the correct distribution. This is obtained with the introduction of an accept/reject test of the Metropolis type, making the procedure a “quasi-heatbath”. The method is described in the next section. Generalizations, including a local version, are presented afterwards. Adler’s stochastic over-relaxation [3] in principle allows to reach z = 1. In practice, this requires careful tuning of the over-relaxation coefficient, depending on the spectrum of A, which is usually difficult to achieve, especially when A itself is fluctuating.